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ABSTRACT 



The time-dependent spectral profile of a resonance line in a homogeneous expanding 
medium is studied by numerically solving an improved Fokker-Planck diffusion equation. 
The solutions are used to determine the time required to reach a quasi- static solution near 
the line center. A simple scaling law for this relaxation time is derived and is fitted to the 
numerical results. The results are applied to the case of Lyman alpha scattering during 
primordial recombination of hydrogen. For a wide range of cosmological models it is found 
that the relaxation times are smaller than the recombination timescale, although not by a 
very large factor. Thus the standard assumption of a quasi-static solution in cosmological 
recombination calculations is reasonably valid, and should not cause substantial errors in 
the solutions. 

1. INTRODUCTION 

The time development of resonance line profiles has evoked interest over the years 
because it plays a role in a variety of astrophysical problems, from the study of the 
resonance excitation of the hydrogen 21 -cm line (Field, 1959) to analyses of the emission 
of x-ray resonance lines (Basko 1978). Previous studies have examined the case of an 
infinite, homogeneous medium that is static, that is, with no macroscopic velocity fields. 
These studies have considered the scattering to be described by a variety of redistribution 
mechanisms (using the notation of Hummer 1962), including Rj redistribution by Field 
(1959), complete redistribution over pure doppler profiles by Ivanov (1967), complete 
redistribution over pure Lorentz profiles by Basko (1978), and type Rjj redistribution 
in the Voigt wings by Basko (1978). These studies of static media have shown that 
the detailed time development is very strongly dependent on the type of redistribution 
mechanism operating in the particular resonance line. 

A non-static case of considerable interest is that of a infinite homogeneous expanding 
medium, since this is a good approximation to the conditions of the universe at the time 
of hydrogen recombination. Provided that the expansion is sufficiently slow, one expects 
that the time-dependent line radiation field will eventually approach an equilibrium (quasi- 
static) configuration. These equilibrium solutions have been extensively studied (see, e.g., 
Rybicki and Hummer 1992 and references therein). In previous studies of the cosmological 
hydrogen recombination epoch (Peebles 1968; Zel'dovich, Kurt and Sunyaev 1969; Krolik 
1989, 1990), it has been assumed that the Lyman alpha line profile rapidly reaches this 
quasi-static configuration, and only the quasi-static solution has then been considered in 
solving the recombination problem. 

The question of whether the Lyman alpha line is in quasi-static equilibrium is very 
important for some aspects of cosmological recombination. The reasons are similar to 
those for the analogous problem of recombination in diffuse nebulae (see, e.g., Osterbrock 
1989), where the trapping of Lyman alpha radiation is an important process in controlling 
the recombination rates, and it is the lack or presence of such trapping that distinguishes 
between Case A and Case B recombination. In the cosmological case, if there is insufficient 
time to build up the quasi-static solution, then the trapping is less effective than otherwise, 
and recombination will be speeded up; this would lead, e.g., to a lower residual ionization 
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of the IGM, which would affect subsequent molecule formation (Palla, Salpeter, & Stahler 
1983; Lepp & Shull 1984). 

The present investigation was motivated by some preliminary, rough estimates of the 
relaxation time for the Lyman alpha line profile, which indicated that it may be comparable 
with the cosmological expansion timescale. If the two timescales were equal, or indeed 
if the relaxation timescale were longer than the expansion timescale, then the previous 
recombination calculations would be in error. The purpose of this paper is to examine in 
detail the relation between these two timescales. 

In §2 the basic equations describing the time dependence of a resonance line profile 
in an expanding universe will be derived. The redistribution mechanism is assumed to be 
Rji, appropriate to the hydrogen Lyman alpha line during the cosmological recombination 
epoch. The Fokker- Planck method is used here for treating Rjj redistribution. The standard 
formulations of this method are due to Harrington (1973) and Basko (1978), following the 
pioneering work of Unno (1955). In this paper we derive and use a new version of the 
Fokker-Planck approximation that has improved behavior near the Doppler core (see the 
Appendix). In §3 analytic results for time dependent development of line radiation fields in 
static media are reviewed, and some new results for the uniform expansion case is given. 
In §4 the numerical method is described and results are given for the relaxation time. A 
simple scaling law is derived that fits the numerical results for the relaxation times quite 
well. These results are then applied to the cosmological recombination problem. Section 
5 contains an evaluation of the importance of some physical effects not treated in our 
equations. Section 6 gives our conclusions. 

2. BASIC EQUATIONS 

For resonance lines with sharp lower levels, such as the Lyman transitions in hydrogen 
at low densities, collisional broadening is negligible and natural broadening dominates. In 
this case the absorption line profile is a Voigt function, but the emission profile depends 
on the spectrum of the absorbed radiation. In the rest frame of the atom, the re-emitted 
(scattered) radiation is coherent, but in the observer's frame the scattering is noncoherent, 
because of the Doppler effect of the atom's thermal motion. This corresponds to type Rjj 
redistribution (Hummer 1962). 

The transfer equation for resonance scattering in an infinite, static, homogeneous and 
isotropic medium can be written as 

l dJ p {vt p ) = _ x<p ^j ^ g +x I z n ^ u <)j [y^ t )du' + C p {t p )ip{u), (1) 
c dtp J 

where J p (u, t p ) is the mean intensity (defined here in terms of photon numbers rather than 
energies) as a function of frequency u and time t p . [The subscript "p" is used here to 
denote "physical" variables that will be replaced by rescaled versions later.] Rn(v,v') 
is the angle- averaged redistribution function describing how the photon absorbed at u' is 
re-emitted at v. (The use of the angle- averaged form is permissible because of the isotropy 
of the radiation field.) Neglecting stimulated emission, the line absorption coefficient is 
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given by XPi"), where 

hi/Q we 2 

X = —n\Bi2 = nifn, (2) 

47r m e c 

and <p{y) is the line profile function (with normalization J <p{y) dv = 1), which is assumed 
to be a Voigt profile. In Eq. (2), e and m e are the electron charge and mass; c is the speed 
of light; n\ is the number density of the lower state; and B\i and /12 are the associated 
Einstein coefficient and oscillator strength. The function C p (t p ) describes the rate of 
creation of new resonance photons due to the radiative cascade following recombination; 
the actual form of C p (t p ) can only be determined on the basis of detailed physical models 
of the recombination process. We have assumed that the newly created photons are emitted 
with the profile <p, which is reasonable for the radiative cascade model. 

Equation (1) applies to a static medium. In an isotropically expanding infinite medium, 
this equation needs to be modified by adding terms to the left hand side corresponding to 
the density change and redshifting of the radiation. Following (Peebles, 1968), these terms 
can be written as 

4^-4^ < 3 > 

where the R(t p ) is the cosmological scale factor; the value of R(t p ) at some standard epoch 
(in our case, the present) will be denoted by R . [Note that the factor 2 (not 3) appears in 
Eq. (3), because of our definition of J p in terms of photon numbers rather than energies.] 
Equation (1) then becomes 

c dt p R dv 

= -XP^Jpfatp) + X J Rn{v,v')Jp(v',tp)dv' + C p {tp)<p{i/). (4) 

It is convenient to change to a dimensionless frequency variable x which measures the 
frequency relative to line center vq in units of the Doppler width, i.e. x = [y — vq) / Az/p 
where Avjj = i> vt /c. Here c is the speed of light and vt is the thermal speed of the atom, 
given by vt = (2ksT /M) 1 / 2 , where ks is Boltzmann's constant, T is the temperature, 
and M is the mass of the atom. For cases of interest here the thermal speeds of the atoms 
are very non-relativistic, vt/c <C 1, so that Eq. (4) becomes 

dJ p (x,t p ) _ cHdJp 

dtp vt ox 

= —kc(j)(x)J p (x,tp) + kc J Rii(x,x')J p (x' ,t p )dx' + Cp(t p )(j)(x), (5) 

where H(t p ) is the local Hubble constant H(t p ) = R(t p ) / R(t p ). The line absorption 
coefficient is now given by k(f>[x), where k = x/Az/£> and (f>(x) is the standard, normalized 
Voigt function, and Ru(x, x') is the standard normalized redistribution function, (see, e.g., 
Rybicki and Hummer 1992, Eqs. [3.1]— [3.5]). Actually we should write (f>(x, a) instead of 
(f>(x) , since the Voigt profile depends on the Voigt parameter a, which is the ratio of Doppler 
to natural width of the line a = A2i/47rAz/£>, where A21 is the Einstein ^-coefficient 
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for the transition. However, for notational simplicity we suppress the Voigt parameter in 
4>(x). For cases of interest here, the Voigt profile is dominated by the Doppler core near 
line center, where <f>(x) « 7r -1 / 2 exp(— x 2 ), so that k is approximately 7T 1 / 2 times the line 
center absorption coefficient. It is convenient to define a new dimensionless time variable 

t = kctp, (6) 

which measures time in units of the mean free time, defined in terms of the integrated line 
absorption coefficient k. We also introduce the rescaled variables 

J{x,t) = {R{t p )/Ro) 2 J p {x,t p ), C{t) = {R{t p )/R ) 2 C p {t p ). (7) 

The transfer equation then takes the form 

^ = -<j>{x)J{x,t) + J ' R n {x,x')J{x',t)dx' + 1 ^ + C{t)cf>{x). (8) 

Here we have introduced the well-known Sobolev parameter 

H _ 8tvH 

1 ~^k~ 3A 2 i A 3 m' ^ 

which characterizes the expansion rate. The second form, in terms of the Einstein A- 
coefficient and wavelength A for the Lyman alpha transition, follows from the Einstein 
relations and the ratio ^2/^1 = 3 for the statistical weights. Physically, 7 is the ratio of 
the mean free path (defined in terms of k) to the distance over which the expansion of the 
medium induces a velocity difference equal to the thermal velocity. Another interpretation 
is that the total optical depth encountered by a photon that redshifts completely through the 
line is ts = sometimes called the Sobolev optical depth. We see that the non-static 
equation is now of the same form as the static equation, except for the additional term 
l(dJ/dx), which describes the redshifting of photons due to the general expansion of the 
medium. 

The treatment of transfer problems involving the redistribution function Ru can 
be quite difficult, since they involve solving integro-differential equations. Fortunately, 
assuming the J{x,t) varies sufficiently slowly over an interval Ax « 1, the scattering 
integral involving R n may be approximated by a differential operator (Harrington 1973; 
Basko 1978, 1981); this is also known as the Fokker-Planck approximation. In the appendix 
we derive an improved version of the approximation, in which the required moments of the 
redistribution function are found exactly. Our approximation is 

/ Rn(x,x')J(x')dx' = cj>(x)J(x) + (^>(x)?^pj , (10) 

which differs from the previous versions in that the full Voigt profile is used in the 
differential term, rather than its asymptotic form 4>(x) ~ a/nx 2 , for \x\ » 1. For this 
reason our approximation is expected to have improved properties near the Doppler core. 
Using Eq. (10), the transfer equation (8) now simplifies to 

dJ(x,t) Id ( .dJ(x,t)\ dJ 

1 ] ~ [ 4>{x) \ ] \ + 1 — + C{t)cj>{x). (11) 



dt 2 dx \ dx dx 
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Apart from the function C(t), this equation depends only on two parameters, the Sobolev 
parameter 7 and the Voigt parameter a. 



It should be noted that the Fokker-Planck approximation to the transfer equation, 
although computationally simple, has some deficiencies. In particular, it is only a good 
approximation when J is not a rapidly varying function of x (since we are neglecting the 
derivatives of order higher than two). The regions in x and t where our solution is not 
reliable will be discussed later. The formulation of Eq. (11) neglects the recoil of the atom 
during scattering, but, as we show in §6, this should not seriously affect our results. 



3. ANALYTIC RESULTS 

There are few known analytic results concerning time-dependent line transfer with 
Rii redistribution. We review and generalize the results of Field (1959) and Basko (1979) 
for static media in sections 3.1 and 3.2. A new exact result for time-dependent coherent 
scattering in an expanding medium is given in §3.3. Finally, in §3.4, we present a simplified 
transfer equation for time-dependent Rjj redistribution, which is applicable for small 7, 
and which has simple scaling properties. 



3.1. Static Medium; The Field Solutions 



Field (1959) found exact solutions for the time development of a radiation field in a 
static medium for Rjj redistribution with a = (also called Rj redistribution), which is 
also a good aproximation when a <C 1 for sufficiently early times. The Field solution 
for an initial field in the form of a Doppler profile J{(x, 0) = 7r -1 / 2 exp(— x 2 ), with zero 
source, C = 0, can be written, after integration by parts and other manipulations, 

In 2 2 /' CO 2 2 

Ji(x,t) =Tv- 1 ' 2 e- x exp(-te~ x )+t e~ w exp(-te~ w )erf (w) dw, (12) 

' J X 

where erf denotes the well known error function (Field 1959; Eq. [16c]). The Field solution 
for a unit, constant source, C(t) = 1, with an initially zero radiation field, is given by 



J s {x,t) 



7T 



-1/2 



00 9 

e w 



exp(-£e x ) 

1 - ( 1 + te~ w " ) exp{-te- w ") erf {w) dw, (13) 



(Field 1959; Eq. [19b]). 



3.2. Static Medium; The Basko Solutions 

Basko (1978) found solutions for Rjj redistribution using a Fokker-Planck approxi- 
mation due to Harrington (1973). The Basko solution for the development from an initial 
pulse of radiation, J(x, 0) = 6(x), and zero source, C(t) = 0, is, 

JAx,t) = - ?—j- e- x4 r\ (14) 
n 1 r(l/4)* 1 /4 ' ^ 1 

where t is a characteristic time t = 8at /it. Although Basko did not give the result J s (x, t) 
for a constant source, C(t) = 1, this can be easily derived from the result (14), if we use 
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the delta-function injection approximation, <f>(x) ~ S(x), for injection of photons in the 
the source term. In this case, J s can be expressed as a superposition of pulses at different 
times, which through a series of straightforward manipulations can be shown to be 

j s (x,t) = J*Ji{x,t?)dt? = 4o r^ 1/4 ) l g l 3r (- 3 / 4 ' z4 /*)- ( 15 ) 

Here r(— 3/4,x 4 /t) denotes an incomplete gamma function (see, e.g., Abramowitz and 
Stegunl964). 

The Fokker-Planck approximation as used by Basko has some limitations: it is valid 
only in the wing regions and it assumes that all the new radiation is injected solely at line 
center, x = 0. These limitations will be clear later when comparing with results from our 
numerical solutions. 

Some of the deficiencies of the Basko approximation could in principle be overcome 
using the analytic results of Grachev (1988), who gave the Green's function for the 
injection of source radiation at arbitrary values of frequency. However, the solution is 
very complicated, involving an integral over Bessel functions. This would have to be 
further integrated over a Voigt profile and integrated over time to obtain the solution for the 
constant source problem. 

3.3. Expanding Medium; Coherent Scattering 

It is possible to give an analytic solution to the non-static equations for the special case 
of no diffusion, that is when the Fokker-Planck term is absent. This does not represent a real 
physical situation, but it provides an interesting and useful reference solution for the more 
realistic cases of R n redistribution. This case is equivalent to that of coherent scattering, 
where we formally substitute R n (x,x') = 4>(x)6(x — x'). Without the diffusion term, 
Eq. (11) becomes, 

dJ dJ „ , 

m=m + c *- (16) 

This equation can easily be solved by changing from variable x to the characteristic 
coordinate x + it. We omit the details, but it is easily checked that the solution is, 

J{x, t) = J{x + 7*, 0) + CV 1 - $(x + it)} , (17) 

where 

[■CO 

= / <t>{x')dx'. (18) 



This function actually depends on the Voigt parameter a (as does (f>(x)), but we again 
suppress this dependence in the notation. The first term of Eq. (17) gives the contribution 
from the initial values of J at t = 0, while the second term gives the contribution from the 
constant source term. 

For the purposes here, we are particularly interested in the case of zero initial field 
J(x, 0) = 0. Then Eq. (17) can be written 

J(x, t) = CV 1 - $(x + it)} . (19) 
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This equation has a simple asymptotic solution for x » 1, where the profile function can 
be approximated by 4>(x) ~ a/nx 2 , so that 



■J(M) ( aC * tV ( 2 °) 

It will be shown later that this provides a good description of the asymptotic region x » 1 
even when diffusion is present. 

3.4. Expanding Medium; Scaling Solution 

It is easy to write an approximate time-dependent transfer equation applicable to 
transfer in the Lorentz wings, in the spirit of Harrington (1973) and Basko (1979). We 
make the replacements 4>(x) ~ a/nx 2 in the diffusion term and use the delta-function 
injection approximation in the source term. Equation (11) then becomes 

dJ(x,t) a d ( 1 dJ\ dJ „ c . . 

The time-independent (quasi- static) solution to this equation is very simple: 

J{x,oo) =C 7 - 1 ( 1 ' ( 3 / o^ l l X 3n (22) 
v ' ' 1 { exp(— 2nqx /3a), if x < 0, v ' 

(Chugai 1980). Note that the solution to the blue (x < 0) is constant, but the solution to the 
red (x > 0) is cut off very sharply at a characteristic frequency of order x c = a 1 / 3 ^ -1 / 3 . 

We are unable to solve the full time-dependent equation (21) analytically; however, it 
is easy to derive some important scaling properties of the solution. Let us introduce the 
new scaled frequency and time variables £ and r defined by 

x = a l l z 1 - ll H, t = a 1 /S" 4 / 3 r. (23) 

It is straightforward to write Eq. (21) in terms of the new variables £ and r. With 
F(£,t) = C~ 1 ~/J(x,t) this equation becomes, 

dF(£,T) 1 d ( 1 dF\ dF r/>N , > 

This equation no longer contains the parameters a, 7 and C, so its solution can be written 
as a function of the scaled variables £ and r alone. These scaling properties may be 
emphasized by writing 

J(x,t;a, 1 ,C) = C 1 - 1 F(Z,T). (25) 

This form will be important for our discussion of relaxation time in the next section. 

Note that Eqs. (21) and (24) only describe frequencies well outside the Doppler 
core. This implies that they are reasonable approximations only when the characteristic 
frequency scale x c = (7/a) 1 / 3 satisfies x c ~> 1. Therefore the condition of validity for 
these equations is 7 <c a. 



4. NUMERICAL RESULTS 



The Fokker-Planck, diffusion-type equation represented by Eq. (11) was solved 
numerically using a straightforward generalization of the the Crank-Nicholson discretization 
scheme (see, e.g., Fletcher 1988; Press et al., 1992, p. 840), modified by the presence of 
the advection term ^/dJ/dx. Being semi-implicit in time, the Crank-Nicholson method has 
good stability properties. Most effort was concentrated on the pure constant source problem, 
C = 1, with zero inital field, J(x,0) = 0, because of its importance to the relaxation 
time problem. The ranges of interest for the parameters 7 and a for the Lyman alpha 
transition during recombination are roughly —10 < log 7 < —6 and —3.5 < log a < —2.5. 
We actually studied a wider range of parameters that included the above ranges, namely, 
— 11 < log 7 < — 1 and —6 < log a < — 1, supplemented by the special cases 7 = (static 
case) and a = (pure Doppler case). 

4.1. Static Medium 

As a test of our method, we applied it first to the static case (7 = 0), for which there are 
analytic solutions. Fig. 1 shows the time-dependent solution of the Fokker-Planck equation 
(11) for the special case of pure Doppler (a = 0), or Rj, redistribution. The exact Field 
solution (13) is also plotted for comparison (dashed lines). The Fokker-Planck equation 
is not expected to be a good description in this case, since the radiation field varies so 
rapidly with x. Indeed the agreement with the Field solution away from line center is 
very poor. However, the agreement near the line center is surprisingly good, since both 
solutions are very flat and of nearly the same level. It is comforting that, even in this most 
inauspicious case, the line center region (which determines the relaxation time) is not too 
badly represented by our improved Fokker-Planck equation. 

Another static (7 = 0) test of the Fokker-Planck equation was a case for which 
a = 10~ 3 , plotted in Fig. 2. At early times, t < 10 4 (not shown), the solutions look 
very much like those of Fig. 1 . At later times the solution begins to broaden considerably 
in the core region and to develop "wings," due to emission in the Lorentz portion of 
the Voigt profile. Eventually, for t » 10 6 , the solutions near line center become very 
well represented by the approximate Basko solution (dashed lines) of Eq. (15). However, 
because of its delta-function injection approximation, the Basko solution does not match 
the solution in the region \x\ » x c , which is dominated by emission into the Lorentz wings. 

4.2. Expanding Medium 

In an expanding medium (7 > 0), the effect of the redshifting of the photons, described 
by the q(dJ/dx) term in Eq. (11), eventually becomes important. The photons that scatter 
significantly outside the doppler core, where the time betweeen scatterings is large, are 
most affected by this redshifting. For those photons scattered blueward of line center, the 
effect of the expansion is to bring them back into the core of the line, where the scattering 
process can begin anew. Red photons, however, are redshifted further away from the line, 
and they are eventually lost. For asymptotically large times this leads to a time-independent 
solution (for a constant source) in which the injection of new photons near line center is 
balanced by the loss due to the redshift (see Rybicki and Hummer 1992). 
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Most of the preceeding effects can be seen in the simple (but otherwise unrealistic) 
case of coherent scattering. In Fig. 3 we plot the coherent solution (19) as log 7 J vs. x for 
a = 10~ 3 and for various values of qt. The line marked "00" is the quasi-static solution 
(~jt — > 00), given by 7 J(x, 00) = Because of the scaling of J and t, Fig. 3 applies to 

all values of 7. It is seen that for early times, 7* <c 1, the solution is approximately equal 
to ~ft(j>(x) and increases linearly with time. As the solution approaches the quasi-static 
form at line center, 7* > 1, the expansion term begins to dominate and the solution skews 
toward the red. Eventually the solution saturates to 7 J ~ 1 in a broad region that extends 
from near line center to a limit in the red which continues to move at a roughly uniform 
rate, much like a wavefront. 

The solution of the scaled Fokker-Planck equation (24) was found numerically using 
the Crank-Nicholson method. The results are plotted in Fig. 4. Due to the suppression 
of the Doppler core in the scaling equation, the early behavior is well described by the 
static Basko solution. At scaled times r ~ 1, the solution begins to skew toward the red, 
and at larger values of r, the solution approaches the time-independent Chugai solution 
(marked "00"). As in the case of the coherent solution, the late-time solution behaves like 
a front moving uniformly towards the red, but now this front is much less sharp. Also 
note the absence of emission in the Lorentz wings, because of the delta-function injection 
approximation implicit in equation (24). 

We next found numerical solutions to the full Rjj redistribution problem, described by 
the Fokker-Planck equation (11). Solutions were obtained for a wide range of parameters, 
but detailed solutions will only be presented for a few representative examples, chosen 
for their relevance to the cosmological recombination problem. In Fig. 5 is plotted the 
numerical solution for a case where a = 10~ 3 and 7 = 10~ 4 . For early times, before the 
effects of expansion become important, the solution looks very much like those in Figs. 1 
and 2, clearly showing the Doppler core. Later the expansion begins to skew the solution 
toward the red and eventually radiation is lost from the line and moves redward in the form 
of a front. Unlike the solution for the scaled equation in Fig. 4, Fig. 5 contains features due 
to emission in the Lorentz wings, and the asymptotic solution is still very much influenced 
by the Doppler core; in particular, the Chugai solution does not describe the asymptotic 
solution very well. The poor description in terms of the scaled Fokker-Planck equation 
(24) is explained by the only marginal satisfaction here of the condition -7 < a. 

In Fig. 6 is plotted a case where again a = 10~ 3 , but now 7 = 10~ 10 , about the most 
extreme case of relevance to the cosmological problem. The frequency spreading is now 
very much larger than before, extending out to hundreds of Doppler widths, and the time 
scales are very much greater. The asymptotic solution is now well represented by the 
Chugai solution (22). In fact, the whole general behavior here is well represented by the 
scaled solution of Fig. 4, except for the very small emission in the wings. 

In the very blue region x » x c , the whole time dependent approach to the asymptotic 
solution is well described by the coherent solution (19), and in particular by the asymptotic 
form (19). This is because diffusion is here completely dominated by the redshifting. 
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5. THE RELAXATION TIME 



5.1 General Theory 



As explained in the Introduction, we are particularly interested in finding out when 
the quasi- static solution to the transfer problem is a good approximation to the full time- 
dependent equations. The quasi-static solution is defined as a time-independent solution 
with the values of the various slowly- varying parameters in the equations, such as C(t), 
^(t), etc., frozen at their local values at some fiducial time. We can investigate this by 
determining the time required for an initially zero radiation field to become approximately 
equal to its asymptotic form in the constant parameter case, where C(t) = C, i(t) = 7, 
etc. This defines a type of relaxation time, which we denote by t r . If this relaxation 
time is much shorter than the characteristic times for the changes in the parameters, i.e., 
C(t) /C'(t), 7(f) /i'(t), etc., then the quasi-static approach is valid. 

Considerable care must be exercised in defining the criteria for the relaxation time. For 
example, from the above numerical solutions in Figs. 3-6 it can be seen that for any time, 
however large, there is always a sufficiently large negative frequency at which the solution 
is far from quasi-static, because there has not been sufficient time for radiation to redshift 
from the line center region. However, we note that the essential region of the quasi-static 
solution for recombination calculations is that near line center, which determines the 
trapping of the resonance radiation. We therefore define our relaxation time t r in terms of 
the time taken to reach the quasi-static solution at x = 0. To be more precise, we choose 
t r to be such that the radiation field at x = reaches some predetermined fraction / of the 
asymptotic value at x = 0. This definition can be formulated concisely through use of the 
residual intensity r(t), defined as 



The choice of the fraction / is somewhat arbitrary. As we shall see, the approach 
to the asymptotic solution is approximately exponential, so one natural choice might be 
/ = 1/e = 0.368. However, we prefer the somewhat more conservative choice / = 0.1, 
so that t r is the time at which the time-dependent solution is within 10% of its quasi-static, 
asymptotic value, or, put another way, the quasi- static solution gives 10% accuracy. The 
exponential behavior then implies that the quasi- static solution gives 1% accuracy at two 
relaxation times, 0.1% at three relaxation times, etc. 

As a simple example, let us evaluate the relaxation time for the case of coherent 
scattering, using the solution (19). Since in this case ^/C~ 1 J(0,oo) = $(0) = 1/2, the 
residual intensity is simply r(t) = 2$(it), and the scaled relaxation time t r is determined 
by the relation / = 2$ (7^). This depends on the Voigt parameter a, but for a < 1 the 
result is essentially the same as that for a = 0, the pure Doppler case, so that the condition 
reduces to the simple one / = erfc(7^ r ). This leads to the result 



r{t) = 1 - J(0,t)/J(0,oo). 
Then the relaxation time t r is determined by the condition 

r(t r ) = f- 



(26) 



(27) 



-1 



(coherent scattering) 



(28) 
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where K c is a constant that depends on the choice for /. For example, our choice / = 0.1 
implies K c = 1.16. 

Now let us derive the relaxation time for the case of Ru redistribution. If we assume that 
the line transfer can be well described by Eq. (21), then a simple derivation of the relaxation 
time follows from the scaling relation, Eq. (25). Using this relation, we may write the quasi- 
static solution as J(x, oo) = C~/~ 1 F(£, oo), so that the ratio of the time-dependent solution 
at to the quasi- static solution at x = takes the form J(0, t)/J(0, oo) = F(0, t)/F(0, oo). 
Therefore the residual intensity is a function of the single variable r, and the relaxation time 
in terms of r is determined by the condition r(r>) = /. In Fig. 7 the residual intensity is 
plotted using the numerical solution of Eq. (24); for r > 1, it is quite accurately represented 
by the exponential law r(r) = 0.920 exp(— 0.663r). For any choice of /, the relaxation 
time in terms of r can be immediately read from this plot, say t> = K. Then from Eq. (23) 
it follows that the scaled relaxation time for Rjj redistribution is 

t r = tfa^V 473 . (29) 
For our choice / = 0.1 we find K = 3.35. 

The formula (29) depends on the correctness of Eq. (21), which was derived on the 
basis of several approximations: that the transfer process is dominated by diffusion in the 
Lorentz wings, and that all new photons are injected exactly at x = 0. Thus the region of 
validity of the formula needs to be checked by direct numerical computations. 

In Fig. 8 we plot typical results the relaxation times found from numerical solution of 
the full Fokker-Planck equation (11), with our preferred choice / = 0.1. This plot is based 
on 126 independent solutions covering the range —6 < log a < —land— 11 < log 7 < — 1. 
There is a region in the upper part of the diagram for which the relaxation times are well 
represented by the asymptotic formula (29) with constant K = 3.35, in agreement with 
that found for the scaled equation. Note that the lower envelope of the scaling region 
corresponds to the condition 7 <c a discussed at the end of §3.4. 

The relaxation times t r in Eqs. (28) and (29) are dimensionless times defined according 
to Eq. (6). The corresponding physical relaxation time t pr for coherent scattering may be 
written 

Ht pr = K c Pt- (coherent scattering) (30) 

Here we have used the definition (9) for 7 and have defined fa = vt/c, the relativistic 
parameter for the thermal motion of the atom, for which typically /?y <C 1. In fact, during 
the recombination epoch, ^«2x 10~ 5 . Since H~ l is the local Hubble time, which is 
a typical expansion time for the medium, Eq. (30) shows that, if coherent scattering were 
an accurate approximation, the relaxation time will always be very much less than the 
expansion time. 

Similarly, the physical relaxation time for Rjj may be written 

Ht pr = Kf3 T (?-J . (31) 

This differs from the result (30) for coherent scattering in having a different proportionality 
constant and, more importantly, in having the extra factor (a/7) 1 / 3 . It is this latter factor 
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that raises the possibility that the relaxation time might be comparable with characteristic 
times for expansion. 



5.2 Application to the Recombination Epoch 



We now wish to compare the relaxation time found above with the typical recombi- 
nation timescales during the cosmological recombination epoch. We used solutions for 
recombination given by Peebles (1968) as well as our own independent solutions to give 
the values of various physical parameters as a function of time (or redshift). 

Let us introduce a number of time scales in addition to the recombination time defined 
above. First, the expansion time t exp is defined by 



where n& is the baryon density. This is simply related to the Hubble time tu = -H" by 
texp = tff/3. More relevant to the problem at hand is the time scale associated with the 
rapid recombination of hydrogen, since this determines the time scale of the function C(t). 
We thus define a recombination time t rec by 



where n e is the electron density. This time is typically 10% of the Hubble time in the 
middle of the recombination phase, roughly where z « 1100. Another time of interest is 
the scattering time t s = 1/kc (cf., Eq. [6]), which represents the free time for a photon 
between scatterings near line center. 

The cosmological model treated here is defined by current values of the density 
parameter 0, the density parameter due to baryons 0&, and the scaled Hubble constant 
h = i?o/(100 km s _1 Mpc -1 ). The cosmological constant is assumed to be zero; in 
any case this should not affect our results. In Fig. 9 we compare various time scales for 
the recombination problem for the set of cosmological parameters = 1,0 = 0.06, and 
h = 0.5. For this case it is seen that the recombination time is always greater than the 
relaxation time, but the ratio between them is as small as ~ 20 for z « 900. 

The case we studied that gives the smallest factor between these two times, while still 
being fairly "reasonable," has O = 0.2, 0& = 0.06, and h = 0.5. The results for this case 
are shown in Fig. 10. Here the ratio between the recombination time and relaxation time is 
~ 5 for a wide range of redshifts, 1000 < z < 1200. 

Because the relaxation time is uniformly smaller than the expansion time for all of 
our cosmological models, we conclude that a quasi-static treatment of the line transfer in 
Lyman alpha is an adequate approximation. The smallest factor between found between 
these times was 5, which, given the exponential approach found in §5.1, suggests that the 
accuracy of the quasi-static solution will be at least as good as one part in 10 5 . 




(32) 




(33) 
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6. OTHER PHYSICAL EFFECTS 



Not included in the above theory is the slight change in frequency during scattering 
due to recoil of the atom. This effect was discussed qualitatively by Adams (1971), and 
a simple, analytic correction to the scattering term was derived by Basko (1981). Making 
Basko's approximations for our scattering term, we obtain the following alteration of 
Eq.(ll), 

dJ 1 d ( ,dJ , \ dJ „. , , 

Here rj is the mean frequency shift in x per scattering due to recoil, given by rj = 
hi>ol '(Mcvt). For the parameters of interest in the hydrogen recombination problem, we 
have r) pd 5 x 10~ 4 . By approximating dJ/dx « J/x, we see that the additional term in 
Eq. (34) only becomes important for frequencies (2??) _1 « 1000, which do not play 
a significant role in determining the relaxation time. At large times, when the radiation 
field is spread to large red frequencies, the redshifting terms in the equation dominate over 
the small recoil effect. Thus, our neglect of this effect is justified. This has been confirmed 
numerically by solving the equations including the additional term for a number of cases 
covering the parameter range of interest. 

A potentially more important effect is that the Lyman alpha scattering process can be 
interrupted, not by collisions (which are negligible for the cosmological problem), but by 
the ambient thermal radiation field (CMB) acting on atoms. During Lyman alpha scattering 
hydrogen atoms alternate between the IS and 2P levels, but exist in the 2P level only 
very briefly (~ 10~ 9 s), almost always returning to the IS level. However, while in the 
2P level, the atom can be radiatively excited to some higher level by the CMB field. The 
subsequent cascade may bring the atom into the metastable 2S level, from which it most 
likely would be excited again by the CMB, but some fraction could undergo two-photon 
decay to IS 1 . The two-photon loss mechanism would tend to reduce the relaxation time, 
since it shortens the scattering process. On the other hand, if the atom is returned to the 2P 
level, so that Lyman alpha scattering continues, it will start by being emitted in the Voigt 
profile, not as a continuation of R n redistribution. This might change the solution by 
putting more radiation further out in the wings of the line, delaying relaxation. Such CMB 
radiative transitions are expected to be quite unlikely for any given scattering event, but 
might nonetheless affect the transfer process because of the large numbers of scatterings 
experienced by a typical Lyman alpha photon. The evaluation of this effect will be left to 
future work. 

7. CONCLUSIONS 

We have investigated the time-dependent radiative transfer in a resonance line with Rjj 
redistribution, which is appropriate to the hydrogen Lyman alpha line during cosmological 
recombination. A new, improved Fokker-Planck approximation was derived, and several 
new analytic results were presented. One principal result is the derivation of the analytic 
formula, Eq. (29), giving the scaled relaxation time for the line radiation field to approach 
its quasi-static form for Rjj redistribution. This formula is written in terms of physical 
variables in Eq. (31). Extensive numerical solutions of the Fokker-Planck transfer equation 
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were performed in order to normalize and check this relaxation formula. 

In order to test if quasi- static models for the line transfer are adequate during 
cosmological recombination, we compared the relaxation time with characteristic time 
scales for recombination, as determined from detailed models of the recombination epoch. 
For all reasonable cosmological parameters it was found that the relaxation time is uniformly 
smaller than the recombination time, but in some extreme cases by a factor of no more 
than 5. We conclude that that reasonably accurate results can be obtained using quasi-static 
solutions as the basis for cosmological recombination calculations. This is quite fortunate, 
since time-dependent calculations involving both radiation and matter would be much more 
complicated. 



In this appendix we derive an improved Fokker-Planck approximation for the Ru 
scattering operator. The derivation is actually somewhat more general than required for 
this paper, in that we treat the case of an angle-dependent radiation field. The isotropic 
case follows as a special case. 

The three-dimensional form of the scattering operator for line scattering with Ru 
redistribution in static media can be written 



where Ru(x, n; x', n') is the angle dependent redistribution function (see, e.g., Hummer 



In the Fokker-Planck method the integral operator defining the scattering term is 
replaced by an second-order differential operator. This method is appropriate when the 
frequency width of the redistribution function is small compared to the scale of frequency 
variation of the intensity. This is precisely the situation for R n redistribution for transfer 
in the Lorentz wings of the line, and we expect that a Fokker-Planck approach will work 
well there. In fact, Fokker-Planck type equations for Ru redistribution were derived and 
used by Unno (1955) and Harrington (1973). However, the equations derived by these 
authors were restricted to the asymptotic wing regime, and become formally singular in the 
core of the line. Furthermore, these equations assumed the angle- averaged form of the Ru 
redistribution function. 

It is the purpose here to provide a derivation of a simple new Fokker-Planck equation 
for Ru redistribution, one based on exact evaluations of the first three moments of the 
complete angle-dependent redistribution function, after Frisch & Bardos (1981). However, 
it is shown that all of these exact moments cannot be used simultaneously, because they are 
inconsistent with photon conservation upon scattering. By introducing a slight alteration of 
the second moment, we obtain a simple Fokker-Planck operator that does conserve photons. 
In its angle- averaged form, this new operator is identical to that used by Harrington in the 
wings of the line, but has the advantage of being nonsingular in the core. 



APPENDIX 

ANGLE DEPENDENT FOKKER-PLANCK OPERATOR FOR R n 





1962). 
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In order to reduce the scattering term to a differential operator we expand the x' 
dependence of I(r, x', n') in a power series about x: 

d 1 

J(r, x', n') = J(r, x, n') + {x 1 - x) — J(r, x, n') + - (x' - x) 2 — /(r, x, n') + • • • . ( A2) 
Then to second order, 

-!-y dO'^ dx'iE77(x,n;x',n')/(r,x',n') = ^ J dfl' [a (x, n, n')J(r, x, n') 

+ Ai (x, n, n') ^/(r, x, n') + -A 2 (x, n, n') /(r, x, n') , {A3) 

where 



A fc (x,n,n') = J dx' (x' - x) k R n (x,n;x',n'). (A4) 
To evaluate the A& we make use of the double Fourier transform representation of Ru, 

„ /"O0 /"O0 . . , , 

R n (T,n;T',ri) = / dx / dxV ra:+ir * # 77 (x, n; x', n'). (A5) 



It has been shown (Rybicki, unpublished; Heinzel 1981) that this can be expressed in the 
closed form 



i2//(r,n; r',n') = g(n ■ n')e 4 r e 4 r e 



i- 2 -A-' 2 -a|r+r'| c -£rrW_ 



(A6) 



Here ^(n • n') is the phase function for the scattering, which can be taken to be isotropic 
(case A) or dipole (case B), 



(A7) 



(A8) 



g A {n ■ n') = 1; g B {n ■ n') = - [l + (n • n') 2 ] . 

The result (6) can be used to evaluate the Fourier transform of A k : 

A k (r,n,ri) = / dx e lTX A k (x, n, n'). 

J — OO 

First of all we write 

R n {T-T',n-T',n')= dxe lTX dx' e lT ( x ~ x >R H {x, n; x', n'), (A9) 

J— CO J — oo 

so that 



A k (t, n, n') = ( — ) R h (t - r', n; r', n') 



r'=0 



(A10) 



Since we need A^ for only the few values k = 0, 1, 2, the easiest way to calculate them is 
by expanding Ru(t — r', n; r', n') to second order in r', 

E 77 (r-r',n;r',n') =ff(n-n')e-i r2 - a l r le^ r '( 1 - nn ')e-^' 2 ( 1 - nn '), 



ff (n • n')e-i r2 - a l r l |l + ^r(l - n • n')r 



/\2 



(l-n-n')- r 2 (l-n.n') 



J2 



(All) 
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Thus we read off the results 



A (r,n,n')=ff(n.n')e-4' 

2i(r,n,n') = -Urg^ • n')(l - n • nV^ 2 "^ 1 , (A12 ) 
2 2 (r,n,n') = ff (n • n') [(1 - n • n') - I r 2 (l - n • n') 2 



Recognizing that exp(— \t 2 — a\r\) is the Fourier transform of the Voigt profile 4>(x) and 
that multiplication by (—ir) is equivalent to d/dx, we have the exact results, 

A (x,n,n') = g(n-n')<j>(x), 

Ai(x,n,n') = -flf(n • n')(l - n • ri)<f>'{x), ^ 13 ) 
A 2 {x,n,n') = g{n ■ n')(l - n • ri)<f>{x) + ^g{n ■ n')(l - n • n')V (*), 

where the primes here denote differentiation with respect to x. 

Although we now have found exact expressions for A^, k = 0,1, 2, it is very important 
not to use the exact expressions for A\ and A% simultaneously. The reason is that the 
differential operator that represents the right hand side of the transfer equation should 
conserve particles in the conservative scattering case; this requires it to be a divergence, 
which in turn requires that 

M = \ A> 2- {A\A) 

We see that the exact values of A\ and Ai do not obey this requirement. In other words, the 
expansion (A2) and (A3) leads to coefficients that do not ensure exact photon conservation 
when truncated to a finite series. 

Experience with radiative transfer problems involving large numbers of scatterings has 
shown that even small errors in photon conservation can lead to very large errors in the 
solution. On the other hand, if exact photon conservation is satisfied, small errors in the 
coefficients lead only to small errors in the solution. Therefore, we shall make Eq. (A 14) 
exact by making appropriate alterations in either A\ or Ai, or both. A careful examination 
of the exact coefficients (A 13) shows that one particularly simple choice stands out, namely 
to adopt the exact expression for A\ , but to approximate Ai by its first term, 

A 2 (x,n,n') = flf(n • n')(l - n • n')<f>{x). (A15) 



We note that this is an excellent approximation in the asymptotic region \x\ » 1, where 
<t>" /<t> ~ 6/x 2 . Near the core, where x < 1, we simply accept the "error" in Ai as being 
unavoidable in any Fokker-Planck treatment that preserves photon conservation. It would 
perhaps be of some interest to try other alterations of A\ and Ai to see whether they 
offer any special advantages. However, the simplicity of the present choice recommends it 
highly, and we shall consider no other here. 
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For coefficients satisfying the relation (A14), the scattering term (A3) can be expressed 
in the form 



^- J dfl' J dx' R n (x,n;x',ri)I(r, x',ri) 



1 

47T 



dfl' 



1 d ( d ^ 

A (x,n,n')/(r,x,n') + -— (A 2 (x,n,ri)—I(r,x,ri) 



(A16) 



Using our adopted values for the scattering term (A 16) can now be displayed explicitly. 
For example, in the isotropic scattering case, 



^- J dfl' J dx' RuA{x,n;x',n')I(r, x',n') 



cf>{x)J(r, x) 



ld_ 
~2dx~ 



(A17) 



and for the dipole case, 
— J dfl' J dx' RiiB{x,n;x',n')I(r,x',n') = 4>{x)- \j{r,x) + riirijKij (r, x)] 



3 d 



8dx I ^^fo. l J ( T,X } + n i n j K ij{ r , x ) - niHi(r,x) - n^n^y^r, x)] 



dx 



where the moments of the radiation field are defined by 

J(r, x) = ~~ J I{r,x,n) dfl 
Hi(r, x) 
Kij(r, x) 



(A18) 



47T 
1 



47T 



I(r, x, n)rii dfl, 
I(r, x, n)nj-ny dfl, 



(A19) 



Lij k (r,x) = J I(r,x,n)ninjn k dfl. 



For notational convenience, an obvious tensor notation (and summation convention) has 
been introduced. 

In many physical situations, e.g., the cosmological case, the radiation field is isotropic, 
and the above formulas simplify greatly, since the odd moments vanish, Hi = and 
Lijk = 0> an d by symmetry the second moment is given by Kij = (l/3)^y J. In this case 
the scattering term has the same form for either the isotropic or dipole phase function, 



J dx' Ru(x, x') J(r, x') = 4>(x)J(r,x) -+ 
It is this form that is used in the body of the paper. 



2dx V K ' dx 



(A20) 
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FIGURE CAPTIONS 

Figure 1. Static solutions (7 = 0) for pure Doppler redistribution (a = 0), found from 
numerical solution of the Fokker- Planck equation (11). The dashed curves are the exact 
analytic results of Field (1959). 

Figure 2. Static solutions (7 = 0) for Rjj redistribution with a = 10~ 3 . The solid curves 
are the numerical results for the Fokker- Planck equation (11) and the dashed curves are 
based on the analytic Basko solution (Eq. [15]). 

Figure 3. Plot of Eq. (19) for coherent scattering with expansion for a = 10~ 3 . The use of 
scaled variables 7 J{x) and make this plot valid for all values of 7. 

Figure 4. Numerical solution of the scaled Fokker- Planck equation (24). Here the scaled 
frequency and time variables are £ = a -1 / 3 ^ 1 / 3 ^ and t = a _1 / 3 7 4 / 3 £. This plot applies 
to any values of a and 7 within the region of validity of the equation (see text). 

Figure 5. Numerical solution of the Fokker-Planck equation (11) for a = 10~ 3 and 
7 = 10~ 4 . The times are given by values of t$ = t/10 5 . The curves below the 10~ 3 curve 
correspond to t$ = 10~ 4 , 10~ 5 , etc. The asymptotic, quasistatic solution is marked "00." 

Figure 6. Same as Fig. 5, except 7 = 10~ 10 and times are given by values of = t/10 12 . 
The curve marked "C" is the Chugai solution (22). 

Figure 7. The residual intensity r(r) vs. t for the scaled Fokker-Planck equation (24). The 
straight dashed line (barely seen in the upper left) is the exponential fit to the lower part of 
the curve. 

Figure 8. Log-log plot of the relaxation time t r vs. 7 for various values of log a. These 
values are based on numerical solution of Eq. (11) with the choice / = 0.1. 

Figure 9. Plot of various time scales (see text) during recombination epoch with 
cosmological parameters = 1, = 0.06, and h = 0.5. 

Figure 10. Same as Fig. 9, except with ft = 0.2 
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